
//=============================================================================================
// Implementation of Exposure Fusion
//
// This work is described in
//   "Exposure Fusion"
//   Tom Mertens, Jan Kautz and Frank Van Reeth
//   In Proceedings of Pacific Graphics 2007
//
// written by Kuen-Shiou Tsai, July 2013
// e-mail: kuenbutsai@gmail.com
//=============================================================================================

#include "expo_fusion.h"

#define NPYRAMID	4
#define NIMAGE		2

//#define ENABLE_LOG

//#define USE_CANNY
//#define USE_LAPLA
#define USE_STDDEV

#define GAMMA0 2.2f
#define GAMMA1 0.45f
#define GAMMA2 1.0f

void onTrackbar(int in, void* object)
{
	expo_fusion* ef = (expo_fusion*) object;

	ef->build_gamma_table();
}

expo_fusion::expo_fusion()
{
	//------------------------------------------------------------------------------
	//	Desc: gamma_i is interger because OpenCV's track bar can only input 
	//		  integer, and we multiply gamma_i by 100 for display 
	//------------------------------------------------------------------------------
	gamma_i.resize(NIMAGE);

	gamma_i.at(0) = (int)(GAMMA0*100.0f);
	gamma_i.at(1) = (int)(GAMMA1*100.0f);

	if(NIMAGE == 3)
		gamma_i.at(2) = (int)(GAMMA2*100.0f);

	build_gamma_table();
}

//------------------------------------------------------------------------------
//	Desc: Interactive gamma table building
//------------------------------------------------------------------------------
void expo_fusion::build_gamma_table()
{
	vector<float>gamma(NIMAGE);

	gamma.at(0) = gamma_i.at(0)/100.0f;
	gamma.at(1) = gamma_i.at(1)/100.0f;

	if(NIMAGE == 3)
		gamma.at(2) = gamma_i.at(2)/100.0f;
	
	if(gamma_table.empty())
		gamma_table.create(NIMAGE, 256, CV_32FC1);

	for(int n=0; n<NIMAGE; n++)
	{
		for(int i=0; i<256; i++)
			gamma_table.at<float>(n, i) = pow(i/256.0f, gamma.at(n));
	}
}

int expo_fusion::process(Mat &src)
{
#ifdef ENABLE_LOG
	clock_t t0 = clock();
#endif

	createTrackbar( "Trackbar0", "recon", &gamma_i.at(0), 512, onTrackbar,(void *)this );
	createTrackbar( "Trackbar1", "recon", &gamma_i.at(1), 512, onTrackbar,(void *)this );
	
	if(NIMAGE == 3)
		createTrackbar( "Trackbar2", "recon", &gamma_i.at(2), 512, onTrackbar,(void *)this );

	vector<Mat> w(NIMAGE);	//weight map
	vector<Mat> g(NIMAGE);	//gamma corrected image

	//------------------------------------------------------------------------------
	//	Desc: synthesis several gamma corrected image
	//------------------------------------------------------------------------------
	gamma_correc(src, g);

	//imshow("gamma0", g.at(0));
	//imshow("gamma1", g.at(1));

#ifdef ENABLE_LOG
	clock_t t1 = clock();
	printf("gamma correct %d\n",t1-t0);
#endif

	//------------------------------------------------------------------------------
	//	Desc: synthesis and normalize weight maps
	//------------------------------------------------------------------------------
	for(int i=0; i<NIMAGE; i++)
		getWeightMap(src, g.at(i), w.at(i));

	normalize_weight(w);

	//imshow("weight0", w.at(0));
	//imshow("weight1", w.at(1));

#ifdef ENABLE_LOG
	clock_t t2 = clock();
	printf("synthesis weight maps %d\n",t2-t1);
#endif

	//------------------------------------------------------------------------------
	//	Desc: build laplacian pyramid for each image
	//------------------------------------------------------------------------------
	vector<vector<Mat>> pyramid_laplacian(NIMAGE);

	for(int i=0; i<NIMAGE; i++)
	{
		pyramid_laplacian.at(i).resize(NPYRAMID);
		build_laplacian_pyramid(g.at(i), pyramid_laplacian.at(i));
	}

#ifdef ENABLE_LOG
	clock_t t3 = clock();
	printf("build laplacian pyramid %d\n",t3-t2);
#endif

	//------------------------------------------------------------------------------
	//	Desc: build gaussian weight pyramid each weight map
	//------------------------------------------------------------------------------
	vector<vector<Mat>> pyramid_gaussian_weight(NIMAGE);
	for(int i=0; i<NIMAGE; i++)
	{
		pyramid_gaussian_weight.at(i).resize(NPYRAMID);
		buildPyramid(w.at(i), pyramid_gaussian_weight.at(i), NPYRAMID-1);
	}

#ifdef ENABLE_LOG
	clock_t t4 = clock();
	printf("build gaussian weight pyramid %d\n",t4-t3);
#endif

	//------------------------------------------------------------------------------
	//	Desc: pyramid blending
	//------------------------------------------------------------------------------
	for(int i=0; i<NPYRAMID; i++)
	{
		for(int k=0; k<NIMAGE; k++)
		{
			//	per element multiply
			pyramid_laplacian.at(k).at(i) = pyramid_laplacian.at(k).at(i).mul(pyramid_gaussian_weight.at(k).at(i));
		}
	}

#ifdef ENABLE_LOG
	clock_t t5 = clock();
	printf("	p-laplacian multiply %d\n",t5-t4);
#endif

	for(int i=0; i<NPYRAMID; i++)
		for(int k=1; k<NIMAGE; k++)
			pyramid_laplacian.at(0).at(i) += pyramid_laplacian.at(k).at(i);

#ifdef ENABLE_LOG
	clock_t t6 = clock();
	printf("	p-laplacian add %d\n",t6-t5);
	printf("pyramid blending %d\n",t6-t4);
#endif

	//------------------------------------------------------------------------------
	//	Desc: reconstruction
	//------------------------------------------------------------------------------	
	reconstruction(pyramid_laplacian.at(0), recon);

#ifdef ENABLE_LOG
	clock_t t7 = clock();
	printf("reconstruction %d\n",t7-t6);
#endif

	//------------------------------------------------------------------------------
	//	Desc: apply
	//------------------------------------------------------------------------------
	Mat yuv;
	cvtColor(src, yuv, CV_BGR2YUV);
	for(int y=0; y<src.rows; y++)
	{
		for(int x=0; x<src.cols; x++)
		{
			if(recon.at<float>(y, x) > 1.0f)
				recon.at<float>(y, x) = 1.0f;
			else if(recon.at<float>(y, x) < 0.0f)
				recon.at<float>(y, x) = 0.0f;

			yuv.at<unsigned char>(y, x*3) = (unsigned char)(recon.at<float>(y, x)*255.0f);
		}
	}
	cvtColor(yuv, recon, CV_YUV2BGR);

#ifdef ENABLE_LOG
	clock_t t8 = clock();
	printf("yuv %d\n",t8-t7);
#endif

	return 0;
}

//------------------------------------------------------------------------------------------
//	Desc: use Gamma correction to synthesis several input images
//------------------------------------------------------------------------------------------
void expo_fusion::gamma_correc(Mat &src, vector<Mat>&dst)
{
	Mat gray;
	cvtColor(src, gray, CV_BGR2GRAY);

	//imshow("gray", gray);

	for(int n=0; n<NIMAGE; n++)
	{
		dst.at(n).create(src.size(), CV_32FC1);

		for(int y=0; y<src.rows; y++)
			for(int x=0; x<src.cols; x++)
				dst.at(n).at<float>(y, x)   = gamma_table.at<float>(n, gray.at<unsigned char>(y, x));
	}
}

//------------------------------------------------------------------------------
//	Desc: Generate a weight map based on several quality measures.
//------------------------------------------------------------------------------
int expo_fusion::getWeightMap(Mat &src, Mat &gc, Mat &weight)
{
#ifdef ENABLE_LOG
	clock_t t0 = clock();
#endif

	Mat hsv;
	cvtColor(src, hsv, CV_BGR2HSV);

#ifdef ENABLE_LOG
	clock_t t1 = clock();
	printf("	CV_BGR2HSV %d\n",t1-t0);
#endif

	//------------------------------------------------------------------------------
	//	Desc: well-exposedness measurement
	//------------------------------------------------------------------------------
	Mat expo(src.size(), CV_32FC1);
	for(int y=0; y<src.rows; y++)
	{
		for(int x=0; x<src.cols; x++)
		{
			float luma = gc.at<float>(y, x);

			if(luma <0.5)
				expo.at<float>(y, x) = 2.0f*luma;
			else
				expo.at<float>(y, x) = -2.0f*luma +2.0f;
		}
	}

	//imshow("expo", expo);	

#ifdef ENABLE_LOG
	clock_t t2 = clock();
	printf("	exposure mesure %d\n",t2-t1);
#endif

	//------------------------------------------------------------------------------
	//	Desc: contrast measurement
	//------------------------------------------------------------------------------
	
	Mat contr;

#if defined (USE_CANNY)
	Mat tmp(gc.size(), CV_8UC1);

	for(int y=0; y<gc.rows; y++)
		for(int x=0; x<gc.cols; x++)
			tmp.at<unsigned char>(y,x) = (unsigned char)(gc.at<float>(y,x)*255.0f);

	Canny(tmp, contr, 64.0, 72.0);

#elif defined (USE_LAPLA)
	Laplacian(gc, contr, CV_32F);
	contr = cv::abs(contr);

#elif defined (USE_STDDEV)
	Mat mu;
	blur(gc, mu, Size(5, 5));

    Mat mu2;
    blur(gc.mul(gc), mu2, Size(5, 5));

	Mat d;
	absdiff(mu2, mu.mul(mu), d);

    sqrt(d, contr);

#endif

	//imshow("contrast", contr);	
	//waitKey(0);

#ifdef ENABLE_LOG
	clock_t t3 = clock();
	printf("	contrast mesure %d\n",t3-t2);
#endif

	//------------------------------------------------------------------------------
	//	Desc: saturation measurement
	//------------------------------------------------------------------------------
	/*Mat satur(src.size(), CV_32FC1);
    for(int y=0; y<src.rows; y++)
    {
        for(int x=0; x<src.cols; x++)
        {
			unsigned char b = src.at<unsigned char>(y,x*3);
			unsigned char g = src.at<unsigned char>(y,x*3+1);
			unsigned char r = src.at<unsigned char>(y,x*3+2);

			float maximum = max(max(b, g), r);
			float minimum = min(min(b, g), r);

			if(maximum == minimum)
				satur.at<float>(y,x)=0;
			else
			{
				float sum = maximum + minimum;
				float diff = maximum - minimum;

				if(sum<=1.0f)
						satur.at<float>(y,x) = diff/sum;
				else
						satur.at<float>(y,x) = diff/(2.0f-sum);
			}
        }
    }*/

#ifdef ENABLE_LOG
	clock_t t4 = clock();
	//printf("	saturation mesure %d\n",t4-t3);
#endif

	//------------------------------------------------------------------------------
	//	Desc: merge
	//------------------------------------------------------------------------------
	weight.create(src.size(), CV_32FC1);

	float satur_weight = 1.0f;
	float contr_weight = 1.0f;
	float expo_weight = 1.0f;

	if((satur_weight == 1.0f) && (contr_weight == 1.0f) && (expo_weight == 1.0f))
	{
		//------------------------------------------------------------------------------------------
		//	Desc: using pow() is time consuming, so avoid call pow() as possible
		//------------------------------------------------------------------------------------------
		for(int y=0; y<src.rows; y++)
			for(int x=0; x<src.cols; x++)
#ifdef USE_CANNY
				weight.at<float>(y, x) = (float)hsv.at<unsigned char>(y, x*3+1) * (float)contr.at<unsigned char>(y, x) * expo.at<float>(y, x);
#else
				weight.at<float>(y, x) = (float)hsv.at<unsigned char>(y, x*3+1) * contr.at<float>(y, x) * expo.at<float>(y, x);
#endif
	}
	else
	{
		for(int y=0; y<src.rows; y++)
		{
			for(int x=0; x<src.cols; x++)
			{
#ifdef USE_CANNY
				weight.at<float>(y, x) = pow((float)hsv.at<unsigned char>(y, x*3+1), satur_weight) *
										 pow((float)contr.at<unsigned char>(y, x), contr_weight) *
										 pow(expo.at<float>(y, x),  expo_weight);
#else
				weight.at<float>(y, x) = pow((float)hsv.at<unsigned char>(y, x*3+1), satur_weight) *
										 pow(contr.at<float>(y, x), contr_weight) *
										 pow(expo.at<float>(y, x),  expo_weight);
#endif
			}
		}
	}

#ifdef ENABLE_LOG
	clock_t t5 = clock();
	printf("	weight merge %d\n",t5-t4);
#endif

	return 0;
}

//------------------------------------------------------------------------------
//	Desc: normalize all weight maps such that weight values sum to one
//		  at each pixel.
//------------------------------------------------------------------------------
void expo_fusion::normalize_weight(vector<Mat> &m)
{
#ifdef ENABLE_LOG
	clock_t t0 = clock();
#endif

	Mat sum = Mat::zeros(m.at(0).size(), m.at(0).type());

	for(int i=0; i<NIMAGE; i++)
		sum += m.at(i);

	for(int i=0; i<NIMAGE; i++)
		m.at(i) /= sum;

	float avg = 1.0f / NIMAGE;

	for(int y=0; y<m.at(0).rows; y++)
	{
		for(int x=0; x<m.at(0).cols; x++)
		{
			if(sum.at<float>(y,x)==0)
			{
				for(int i=0; i<NIMAGE; i++)
					m.at(i).at<float>(y,x) = avg;
			}
		}
	}

#ifdef ENABLE_LOG
	clock_t t1 = clock();
	printf("	---------\n");
	printf("	normalize %d\n",t1-t0);
#endif
}

//------------------------------------------------------------------------------------------
//	Desc: Lossless Image Decomposition, please refer to "The Laplacian Pyramid as a Compact 
//		  Image Code" Burt, P., and Adelson, E. H., IEEE Transactions on Communication, 
//		  COM-31:532-540 (1983)."
//------------------------------------------------------------------------------------------
void expo_fusion::build_laplacian_pyramid(Mat &src, vector<Mat> &pyramid_laplacian)
{
#ifdef ENABLE_LOG
	clock_t t0 = clock();
#endif

	vector<Mat> pyramid_gaussian;
	buildPyramid(src, pyramid_gaussian, NPYRAMID-1);

#ifdef ENABLE_LOG
	clock_t t1 = clock();
	printf("	buildPyramid %d\n",t1-t0);
#endif

	for(int i=0; i<NPYRAMID-1; i++)
	{
		Mat pUp;
		pyrUp(pyramid_gaussian.at(i+1), pUp);

		pyramid_laplacian.at(i) = pyramid_gaussian.at(i)- pUp;
	}

#ifdef ENABLE_LOG
	clock_t t2 = clock();
	printf("	pyrUp %d\n",t2-t1);
#endif

	pyramid_laplacian.back() = pyramid_gaussian.back().clone();

#ifdef ENABLE_LOG
	clock_t t3 = clock();
	printf("	clone %d\n",t3-t2);
#endif
}

//------------------------------------------------------------------------------------------
//	Desc: Lossless Image reconstruction, please refer to "The Laplacian Pyramid as a Compact 
//		  Image Code" Burt, P., and Adelson, E. H., IEEE Transactions on Communication, 
//		  COM-31:532-540 (1983)."
//------------------------------------------------------------------------------------------
void expo_fusion::reconstruction(vector<Mat> &pyramid_laplacian, Mat &recon)
{
#ifdef ENABLE_LOG
	clock_t t0 = clock();
#endif

	pyrUp(pyramid_laplacian.at(NPYRAMID-1), recon);

#ifdef ENABLE_LOG
	clock_t t1 = clock();
	printf("	first pyrUp %d\n",t1-t0);
#endif

	for(int i=NPYRAMID-2; i>=1; i--)
	{
		recon += pyramid_laplacian.at(i);
		pyrUp(recon, recon);
	}

#ifdef ENABLE_LOG
	clock_t t2 = clock();
	printf("	pyrUp %d\n",t2-t1);
#endif

	recon += pyramid_laplacian.at(0);

#ifdef ENABLE_LOG
	clock_t t3 = clock();
	printf("	last stage %d\n",t3-t2);
#endif
}

//======================================================================================

/*
int expo_fusion::cl_init()
{
	cl_int err = CL_TRUE;

	err = clGetPlatformIDs(1, &platform, NULL);
	if(CL_SUCCESS != err)
	{
		printf("clGetPlatformIDs failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clGetPlatformIDs succeed: platformID: %x found\n", (unsigned int)platform);
	//-----------------------------------------------------------------------------------

	err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
	if(CL_SUCCESS != err)
	{
		printf("clGetDeviceIDs failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clGetDeviceIDs succeed: deviceID: %x found\n", (unsigned int)device);
	//-----------------------------------------------------------------------------------

	ctx = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
	if(CL_SUCCESS != err)
	{
		printf("clCreateContext failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clCreateContext succeed: contextID: %x\n", (unsigned int)ctx);
	//-----------------------------------------------------------------------------------

	cmdq = clCreateCommandQueue(ctx, device, 0, &err);
	if(CL_SUCCESS != err)
	{
		printf("clCreateCommandQueue failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clCreateCommandQueue succeed: commandqueueID: %x\n", (unsigned int)cmdq);
	//===================================================================================

	char* kernelSource = cl_ReadSources("kernel.cl");
	program = clCreateProgramWithSource(ctx, 1, (const char**)&kernelSource, NULL, &err);
	if(CL_SUCCESS != err)
	{
		printf("clCreateProgramWithSource failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clCreateProgramWithSource succeed: programID: %x\n", (unsigned int)program);
	//-----------------------------------------------------------------------------------

	err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
	if(CL_SUCCESS != err)
	{
		printf("clBuildProgram failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clBuildProgram succeed\n");
	//-----------------------------------------------------------------------------------

	kernel = clCreateKernel(program, "gamma_correc", &err);
	if(CL_SUCCESS != err)
	{
		printf("clCreateKernel (gamma_correc) failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clCreateKernel (gamma_correc) succeed: kernelID: %x\n", (unsigned int)kernel);

	return 0;
}

int expo_fusion::cl_initBuffer(Mat &src, Mat &dst)
{
	dst.create(src.size(), CV_32FC3);

	cl_int err = CL_TRUE;

	input_buf = clCreateBuffer(ctx, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(unsigned char)*src.rows*src.cols*src.channels(), src.data, &err);
	if(CL_SUCCESS != err)
	{
		printf("clCreateBuffer [0] failed: error: %d\n", err);
		exit(err);
	}
	else
		printf("clCreateBuffer [0] succeed: bufferID: %x\n", (unsigned int)input_buf);
	//-----------------------------------------------------------------------------------

	output_buf = clCreateBuffer(ctx, CL_MEM_WRITE_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float)*dst.rows*dst.cols*dst.channels(), dst.data, &err);
	if(CL_SUCCESS != err)
	{
		printf("clCreateBuffer [1] failed: error: %d\n", err);
		exit(err);
	}
	else
	{
		printf("clCreateBuffer [1] succeed: bufferID: %x\n", (unsigned int)output_buf);
	}

	return 0;
}

int expo_fusion::cl_execute_kernel(Mat &src, Mat &dst, float gamma)
{
	cl_int err = CL_TRUE;

	err  =  clSetKernelArg(kernel, 0, sizeof(cl_mem), (void *) &input_buf);
    err  |= clSetKernelArg(kernel, 1, sizeof(cl_mem), (void *) &output_buf);
	err  |= clSetKernelArg(kernel, 2, sizeof(cl_float), &gamma) ;

    if (err != CL_SUCCESS)
    {
        printf("ERROR: Failed to set input kernel arguments (%d)...\n", err);
		exit(err);
        return false;
    }

	size_t global_work_size[2] = { src.cols, src.rows};
    size_t local_work_size[2]= { 16, 16};							

    // execute kernel
    if (CL_SUCCESS != (err = clEnqueueNDRangeKernel(cmdq, kernel, 2, NULL, global_work_size, local_work_size, 0, NULL, NULL)))
    {
        printf("ERROR: Failed to run kernel...\n");
        return false;
    }
    
	err = clFinish(cmdq);

	void* tmp_ptr = NULL;
    tmp_ptr = clEnqueueMapBuffer(cmdq, output_buf, true, CL_MAP_READ, 0, sizeof(float)*dst.rows*dst.cols*dst.channels(), 0, NULL, NULL, NULL);
	if(tmp_ptr != dst.data)
	{
		printf("ERROR: clEnqueueMapBuffer failed to return original pointer\n");
		return false;
	}


    err = clFinish(cmdq);

    clEnqueueUnmapMemObject(cmdq, output_buf, tmp_ptr, 0, NULL, NULL);

	return 0;
}

int expo_fusion::cl_release()
{
	clReleaseMemObject(input_buf);
	input_buf = NULL;
	clReleaseMemObject(output_buf);
	output_buf = NULL;

	clReleaseKernel(kernel);
	clReleaseProgram(program);
	clReleaseCommandQueue(cmdq);
	clReleaseContext(ctx);

	return 0;
}

char* expo_fusion::cl_ReadSources(const char *fileName)
{
	FILE *file = fopen(fileName, "rb");
	if (!file) {
		printf("ERROR: Failed to open file '%s'\n", fileName);
		return NULL;
	}

	if (fseek(file, 0, SEEK_END)) {
		printf("ERROR: Failed to seek file '%s'\n", fileName);
		fclose(file);
		return NULL;
	}

	long size = ftell(file);
	if (size == 0) {
		printf("ERROR: Failed to check position on file '%s'\n", fileName);
		fclose(file);
		return NULL;
	}

	rewind(file);

	char *src = (char *) malloc(sizeof(char) * size + 1);
	if (!src) {
		printf("ERROR: Failed to allocate memory for file '%s'\n", fileName);
		fclose(file);
		return NULL;
	}

	printf("Reading file '%s' (size %ld bytes)\n", fileName, size);
	size_t res = fread(src, 1, sizeof(char) * size, file);
	if (res != sizeof(char) * size) {
		printf("ERROR: Failed to read file '%s'\n", fileName);
		fclose(file);
		free(src);
		return NULL;
	}

	src[size] = '\0';
	fclose(file);

	return src;
}

*/

